Neel to dimer transition in spin-S* antiferromagnets: Comparing bond operator theory 
with quantum Monte Carlo simulations for bilayer Heisenberg models 
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We study the Neel to dimer transition driven by interlayer exchange coupling in spin-S Heisenberg 
antiferromagnets on bilayer square and honeycomb lattices for 5=1/2, 1, 3/2. Using exact stochastic 
series expansion quantum Monte Carlo (QMC) calculations, we find that the critical value of the 
interlayer coupling, Jxc[S], increases with increasing S, with clear evidence that the transition is in 
the 0(3) universality class for all S. Using bond operator mean field theory restricted to singlet and 
triplet states, we find JxcfS] oc S{S + 1), in qualitative accord with QMC, but the resulting JxcfS] 
is significantly smaller than the QMC value. For 5=1/2, incorporating triplet-triplet interactions 
within a variational approach yields a critical interlayer coupling which agrees well with QMC. For 
higher spin, we argue that it is crucial to account for the high energy quintet modes, and show that 
including these within a perturbative scheme leads to reasonable agreement with QMC results for 
S'=l,3/2. We discuss the broad implications of our results for systems such as the triangular lattice 
S=l dimer compound Ba3Mn208 and the S = 3/2 bilayer honeycomb material Bi3Mn40i2(N03). 



I. INTRODUCTION 



Spin dimer compounds provide the simplest realiza- 
tion of a magnetically disordered ground state — one 
where strongly coupled pairs of spins entangle to form 
singlets. Such systems are also of great interest since they 
undergo magnetic field induced spin ordering via a quan- 
tum phase transition which is analogous to Bose-Einstein 
condensationFi"— There are many well-known spin dimer 
compounds^"— and well studied model HamiltoniansSrJ^ 
exhibiting such physics for 5*= 1/2 spins. However, on- 
going experiments on higher spin systems, such as the 
S=l triangular lattice dimer compoundl^ Ba3Mn208, 
point to a need to better understand higher spin gener- 
alizations and instabilities of such dimer states driven by 
inter-dinier interactions. 

Here, we explore this issue using a simple model which 
exhibits such a dimerized ground state - the bilayer 
Heisenberg antiferromagnet with a Hamiltonian given by 





H 



J± ^ Si,i • Si_2 + Ji ^ ^ ^i,e ■ Sj 



(1) 



(ij) 1= 



Here, i labels sites in one layer, £ = 1, 2 is the layer index, 
and {ij) represents nearest neighbor pairs of spins within 
each layer. For J_l ^ Ji, the first term in H dominates 
and the ground state is composed of isolated interlayer 
singlets with S^^i -I- Si, 2 = for every i. If J± <^ Ji, the 
system will order magnetically provided the second (in- 
tralayer) term in the Hamiltonian is not too frustrated by 
the lattice geometry. Here, we restrict our attention to 
cases where each layer is itself a bipartite lattice so that 
the ground state for J± ^ Ji has long-range Neel order. 
This model Hamiltonian has been extensively studied for 
the S" = 1/2 square lattice bilayer— li^"—. Effects of dis- 
order, induced by site dilution, have also been exploredJ^ 
However, there has been relatively little work on under- 



1 



FIG. 1: The interlayer dimer state on square and honeycomb 
bilayers with singlet correlations between layers. The in-plane 
exchange Ji and the interplane exchange J± act as shown. 



standing the higher spin generalizations of the Hamil- 
tonian in Eq. [T] The spin-S square lattice bilayer has 
been studied using Schwinger boson mean field theorjj2£ 
and series expansions. Variants of this spin-S model 
have been argued to support novel spin solid phases in 
three dimensionSi-^ In this paper, we study this Hamil- 
tonian for S = 1/2,1,3/2 spins on square and honey- 
comb bilayers using exact quantum Monte Carlo simula- 
tion algorithms'^ and approximate analyses based on a 
bond operator method generalized to arbitrary spini^ 

Our main results are as follows. (i) Using exact 
stochastic series expansion quantum Monte Carlo (QMC) 
calculations, we find a Neel to dimer transition with in- 
creasing J± that is in the 0(3) universality class for all 
the models we have studied, (ii) The critical value of 
the interlayer coupling, Jicl-S'], for the Neel to dimer 
transition is found to increase for higher spin. Using 
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a bond operator mean field theory restricted to singlet 
and triplet states, we find J±c[S] oc S{S + 1), in qual- 
itative accord with QMC results. However, there is a 
quantitative discrepancy between the mean field J_lc[5'] 
and its QMC value, which becomes more significant for 
higher spin, (iii) For S = 1/2, we show that taking 
into account triplet-triplet interactions within a vari- 
ational approach brings the J±c[S] value close to the 
QMC result. For higher spin, we show that the domi- 
nant corrections to the critical point arise from the high 
energy quintet modes and direct triplet-triplet interac- 
tions are less important. Incorporating the quintet ex- 
citations within a perturbative treatment is shown to 
yield a critical interlayer coupling which is in good agree- 
ment with QMC results for S = 1,3/2. We discuss the 
broad implications of our results for high spin antifer- 
romagnets such as the triangular lattice S — 1 dimer 
compound^ Ba3Mn2 08 and the 5 = 3/2 bilayer honey- 
comb materia^^ Bi3Mn40i2(N03). 

This paper is organized as follows. Section II contains 
results from the QMC simulations on the phase diagram 
of the honeycomb and square lattice bilayer models for 
S — 1/2, 1, 3/2. In Section III, we outline the bond oper- 
ator formalism generalized to the case of spin-S. Section 
IV gives bond operator mean field theory results for the 
square and honeycomb lattice models. Section V dis- 
cusses the variational approach that we use to take into 
account corrections beyond mean field theory. Section VI 
analyses the 5* = 1/2 model including the effect of triplet- 
triplet interactions, while Section VII contains a treat- 
ment of the dominant quintet corrections for S > 1/2. 
We end with a discussion in Section VIII. Details are 
contained in Appendices. 



II. QUANTUM MONTE CARLO SIMULATIONS 

The bilayer honeycomb and square lattices are bipar- 
tite lattices which can be split into two sublattices A and 
B with every lattice bond being a link between sites be- 
longing to different sublattices. This ensures that there is 
no sign problem, so that the model in Eq. [T] is amenable 
to quantum Monte Carlo simulations. We perform quan- 
tum Monte Carlo simulations for S — 1/2, 1,3/2 on the 
bilayer square (of linear system size L = 12, 16, 24, 32, 40) 
and bilayer honeycomb {L = 12,18,24,30,36) lattices 
using the Stochastic Series Expansion algorithm;22, For 
S = 1 and S = 3/2, simulations are performed with mod- 
ified worm weights, which lead to a slightly more efficient 
algorithm, as in Ref. IH. At large enough ratio J±/Ji, 
the system undergoes a quantum phase transition from a 
Neel state to a dimerized paramagnetic state. To locate 
quantum critical points, we perform finite size scaling 
analysis of the superfluid density and the staggered mag- 
netization density squared. We measure the superfluid 
density ps by measuring winding number fluctuations^ 




FIG. 2: Scaling of the superfluid density (upper panel) and 
of the staggered magnetization density squared (lower panel) 
for the S — 1 antiferromagnet on the bilayer square lattice. 
The curves cross at a distinct point around Jx/Ji — 7.15. 
The insets show the corresponding data collapse for z = 1, 
V = 0.7112, /3 = 0.3689, and Jxc/Ji = 7.15. Lines guide the 
eye. The error bars are smaller than the symbol size if not 
visible. 



where W\^2 are the winding numbers in two spatial di- 
rections and T is the temperature. The staggered mag- 
netization density squared is given by 



1 ^ 



where (—1)'' = 1 for lattice sites from sublattice A, 
(— 1)P = —1 for sites from sublattice B, and N is the 
number of lattice sites. In the vicinity of a continuous 
phase transition the superfluid density scales as 



= L 



2-d-;. 



F{[K-K,]L''\—) 



T 



2Ji 



where L is the linear system size, d = 2 is the dimension- 
ality of the system, T is the temperature, [K — Kc] = 



3 
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FIG. 3: Scaling of the superfluid density for the S — 3/2 an- 
tiferromagnet on the bilayer square lattice. The curves cross 
at a distinct point around J±/Ji = 13.634. The inset shows 
the corresponding data collapse for z = 1, v = 0.7112, and 
J±c/Ji ~ 13.634. Lines guide the eye. The error bars are 
smaller than the symbol size if not visible. 

[(Jj_ — J±c)/Ji] is the distance from the critical point 
J±c/ Ji, and v is the correlation length critical exponent. 
The staggered magnetization density squared scales as 

where is the critical exponent. If one plots PsL^ as 
a function of J±/Ji at large enough and fixed value of 
\/{TL^) then the curves for different system sizes should 
cross at the critical point Ja.c/Ji- If one plots psL^ as 
a function of [K — K^L^^"^ , with appropriately chosen 
values of the critical exponents and Kc, the curves for 
different systems sizes should collapse onto the universal 
curve given by the function F. Similarly, jm^pL^''/'' as 
a function of J±/ Ji should have a distinct crossing point 
at the critical point and \ms\'^ L'^^^'^ as a function of [K — 
Kc\L^^^ should collapse onto the universal curve given by 
the function M . We perform simulations at fixed aspect 
ratio T = Ji/2L. 



A. Square lattice 

The quantum critical point for the S — 1/2 bilayer 
quantum antiferromagnet on the square lattice was found 
in Refs. [nllli J^c/Ji = 2.5220(1). In the present work, 
we find that the quantum critical points are located at 
Jj^^/ Ji = 7.150(2) for S* = 1, and at J±c/ Ji = 13.634(3) 
for S = 3/2. The data scale very well with the critical 
exponents v — 0.7112 and /3 = 0.3689 of the 0(3) univer- 
sality class— for any value of spin. The crossing points 
and data collapse for 5 = 1 and 5 = 3/2 are shown in 
Figs. [5] and |31 Note that we do not show the scaling of 




4.74 4.76 4.78 4.8 4.82 4.84 




4.74 4.76 4.78 4.8 4.82 4.84 



FIG. 4: Scaling of the superfluid density (upper panel) and of 
the staggered magnetization density squared (lower panel) for 
the S = 1 antiferromagnet on the bilayer honeycomb lattice. 
The curves cross at a distinct point around Jx/Ji = 4.785. 
The insets show the corresponding data collapse for z = 1, 
u = 0.7112, P = 0.3689, and J_ic/Ji = 4.785. Lines guide the 
eye. The error bars are smaller than the symbol size if not 
visible. 



the magnetization density squared for 5 = 3/2 because 
the data points are too noisy. 



B. Honeycomb lattice 

We find that that for the honeycomb lattice the quan- 
tum critical points are located at J±c/ Ji = 1.645(1) for 
5 = 1/2, J^rJJi = 4.785(1) for 5=1, and Jj.c/Ji = 
9.194(3) for 5 = 3/2. The data scale very well with the 
critical exponents v = 0.7112 and /? = 0.3689 of the 0(3) 
universality class^ for any value of spin. The crossing 
points and data collapse for 5 = 1 and 5 = 3/2 are 
shown in Figs. [3] and [HI We do not show the scaling of 
the magnetization density squared for 5 = 3/2 because 
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FIG. 5: Scaling of the superfiuid density for the S — 3/2 
antiferromagnet on the bilayer honeycomb lattice. The curves 
cross at a distinct point around Jx/Ji — 9.194. The inset 
shows the corresponding data collapse for z — 1, v — 0.7112, 
and J±c/Ji = 9.194. Lines guide the eye. The error bars are 
smaller than the symbol size if not visible. 



the data points are too noisy. 
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FIG. 6: J±c/Ji[S] as a function of S{S + 1) for the bilayer 
square and honeycomb lattices. Lines are linear fits. Note 
that the curves cross approximately at S{S + 1) = 0. 



C. Critical point as a function of spin 

In Fig. ini we show the quantum critical points Jj_c/ Ji 
as functions of S{S + 1) for the bilayer square and honey- 
comb lattices. We find that J±c/Ji[S] is a linear function 
of S{S + 1). In the following sections, we will attempt to 
make sense of these QMC results using an extension of 
the bond operator theory of the dimerized state and its 
instability to Neel order. 



III. BOND OPERATOR REPRESENTATION 

An elegant approach which allows us to understand the 
physics of the dimer ground state and its magnetic or- 
dering instabilities is the bond operator formalism which 
was first proposed for 5 = 1/2 antiferromagnetsi^ In 
this scheme, the spin operators are represented in a new 
basis consisting of singlet and triplet states on the inter- 
layer bonds (i, !)-(«, 2). In the limit where the intralayer 
coupling Ji — 0, the ground state consists of localized 
singlets on these bonds, with a gap J± to the triplet 
excitations. A nonzero Ji allows a pair of neighbor- 
ing bonds (i,l)-(i,2) and (j, l)-(j, 2) to exchange their 
singlet/triplet character. Such a 'triplet hopping' pro- 
cess converts the localized triplet modes into dispersing 
'triplons', with three-fold degenerate bands due to the 
underlying SU (2) symmetry of the Hamiltonian. In this 
picture, the dimer to Neel transition is an 0(3) transi- 
tion driven by the condensation of triplon modes at a 
certain wavevector where the dispersion minimum hits 
zero. Generalizations of this approach to spin-1 magnets 
have been proposed earlier .^"'^'^^ Here, we adopt a recent 
generalization of the bond operator method to arbitrary 
spin^i to study bilayer Heisenberg antiferromagnets. 

In a spin-S* bilayer system, in the limit J± 3> Ji, we 
have isolated interlayer bonds. The bond can be in one of 
the following states: a singlet, a 3-fold degenerate triplet, 
a 5-fold quintet, etc. We introduce one boson for each of 
these states: 

^ 4|0), 
|ij,me{-l,0,l}) = 4,m|0)> 
h,me{-2,- .2}) = qLJO), 

The index i here runs over all interlayer bonds, and m 
labels the i'j-component of the total spin on the inter- 
layer bond. These boson operators form the basis for a 
bond operator representation. To restrict to the physical 
Hilbert space of spins, every interlayer bond should have 
exactly one boson, 

4^*+ E 4,,A™+ 9l,„*,n + • • • = 1. (2) 

T?x=-l,0,l ri=-2,---,2 

In terms of bond operators, the exchange interaction on 
an interlayer bond is given by 

m=-l,0,l 
m=~2,--- ,2 

where e, = -J^S{S + 1), = Jj.{l - S{S + 1)}, and 
Eg = J_i{3- 5(5 + 1)}. 

Bond operator theory re-expresses the spin operators 
and their interactions in terms of these bond bosons. In 



the hmit J± ^ Ji, the singlets, triplets, quintets, etc. 
form a hierarchy with the energy spacing between each 
tier of order J±. In this paper, we restrict our analysis to 
the low energy subspace of singlets, triplets and quintets 
on a bond, and neglect higher spin states as they are 
much higher in energy. 

We first turn to the usual bond operator mean field 
theory retaining only singlet and triplet modes, ignoring 
triplet interactions and higher excited states and impos- 
ing the constraint in Eq. [2] on average. We then con- 
sider, in turn, the effect of triplet-triplet interactions for 
S = 1/2 and the effect of quintet states for S > 1/2. For 
convenience of notation, we henceforth set Ji = 1, thus 
measuring J± in units of Ji. 



IV. SINGLET-TRIPLET MEAN FIELD THEORY 

At mean field level, the interlayer dimer state is de- 
scribed by a uniform condensate of the singlet bosons. 



with (sj) 



J\ - 



s. Retaining only triplet excitations. 



the spin operators at each site are given by^l 

+ ^{4,1*^,0 + (4) 

+ \{tllt^^-tl-lt^,-l}■ (5) 

Using these expressions, the Hamiltonian takes the 
form 

i,m i \ m / 



+ 



25(5 + 1) 



E [{t^,o + tlo}{tj,o + tlo} 



{h-l-tl,}{t].-tj^i} + h. 



(6) 



where is a Lagrange multiplier which enforces the con- 
straint in Eq. [5] on average. N± is the number of in- 
terlayer bonds. We have dropped quartic terms in the 
triplet operators (which corresponds to ignoring triplet- 
triplet interactions). 

In the rest of this paper, we use the following two ba- 
sis sets to represent triplet states: \to)i, \ti)i} or 
{\tx)i, \ty)i, \tz)i}- The former basis labels states by the 
z-projection of spin. The latter labels each state by the 
direction in which its spin projection is zero. We can 
go from one basis to another using \to)i = \tz)i and 
\t±i}i = {T\tx}i — 'i\ty)i)/V2- Below, we will use the 
index m to represent an element of the first basis and u 
to represent an element of the second. 




FIG. 7: Top view of bilayers. (Top) Square lattice with prim- 
itive lattice vectors x and y shown. (Bottom) Honeycomb lat- 
tice. The shaded region is the unit cell composed of two sites. 
Sites marked with a red circle belong to the A sublattice. Un- 
marked sites belong to the B sublattice. The primitive lattice 
vectors a and b are shown. 



A. Square lattice bilayer 

A top view of the square lattice bilayer with the rele- 
vant primitive lattice vectors is shown in Fig.[7l At mean 
field level, the Hamiltonian of Eq. [T] may be written as 



=-J±N^S{S + l)s^ 



E' 

lii.u^{x.y.z} 



- ns^N± - 
2ek 



2ek 



3N,A 



^k,«,(7) 



where V'k.u = [^k,ti ^-kul"^- primed summation 

indicates that if k is included in the sum, then — k is 
excluded. The coefficients in the Hamiltonian matrix are 



A = Ji{l-5(5+l)}-/i, (8) 
2S(S + 1) 

Ck = — 4 ^s^(cos(fc^) +cos(fcy)). (9) 



Diagonalizing this Hamiltonian matrix by a bosonic Bo- 
goliubov transformation (see Appendix we obtain 
eigenvalues Ak = ^/A(A'^\-4£^ for the energies of the in- 
dependent 'triplon' modes. Each of these modes adds a 
zero point contribution to the ground state energy, yield- 
ing 



3iVj.yl 



3E'A. 



(10) 
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We minimize this ground state energy with respect to fi 
and s, via dE^^^/d^ = and dE^°^/ds'^ = 0. This yields 
the two equations 



Aw ' 



-J^SiS + l) 



-T 



(11) 
(12) 



Using the values of s and /i thus obtained, we may calcu- 
late the gap to triplet excitations. The dimer-Neel tran- 
sition occurs when the triplon gap vanishes at J± = J±c- 
We have explicitly checked that triplon condensation 
at k = (tt, tt) yields Neel order on the bilayer. Using 
Egns. fTlll2l above, we arrive at the following two results 
at the critical point, (i) The value s at the dimer-Neel 
critical point is independent of spin and is given by 



2Ni 



4 + (cos fca: + cos ky ) 

± y/A + 2(cosfca; -I- COSfcj^) 



E 

k 



(13) 



A numerical evaluation shows Sc ~ 0.904. (ii) We find 
the location of the dimer-Neel critical point 



40 32 ^ , 



1 



' 2 (cos kx + cos ky] 



(14) 



A numerical evaluation yields J±c ~ 3.047S'(S' + 1). For 
5 = 1/2, this mean field result, J±c[S = 1/2] w 2.286, 
agrees with previous workii and is slightly smaller than 
the QMC valueJi For higher spin, the mean field esti- 
mates, JjLc[S = 1] « 6.095 and Jjlc[S = 3/2] « 11.428, 
are significantly smaller than our QMC results. This 
comparison is summarized in Table ID The scaling result 
J_Lc ~ 5'(5' -1-1) has been suggested in Ref. ,21 on the basis 
of series expansion calculations. Remarkably, as shown 
in FiglHl this scaling relation derived from mean-field the- 
ory seems to be reasonably accurate even for exact QMC 
results. 



B. Honeycomb lattice bilayer 

The honeycomb lattice is composed of two interpene- 
trating triangular lattices, as shown in Figl?) Operators 
therefore come with an additional sublattice index which 
distinguishes A and B sublattices. The mean field Hamil- 
tonian is given by 



rr(0) _ 

- 



3N^C 



(15) 



where C = (J±{1 - S{S + 1)} - fi). N±_ denotes the 
number of interlayer bonds in the honeycomb bilayer. 
The operator ipk.u and the Hamiltonian matrix Mk are 



given by 



k, (i 



/ ik,A,u \ 

t-k.B,u 

t 

-k,A,ti 
V ^-k,_B,M / 



A/k = 



C /3k /3k 

/3k C /3* 

/3k C /3k I '^'"'^ 

/3* /3* C 



^ — ikb 



-, — ika—if^b 



where ^k = 2^^^^s^-fk, with 7k = 1- 
and we have defined ka = k • d and fcf, = k • &. Diago- 
nalizing this Hamiltonian (see Appendix |B]) , we obtain 
two eigenvalues for every k. The eigenvalues are given 
by Ak.i/2 — \/C^ =F 2C|/3k|. The mean field ground state 
energy is given by 



E. 



(0) 

O 



3Nj_C 



k,l + Ak,2) 



(17) 



As before, we demand dEjj /d^ ~ and dE^^ /ds^ — 0. 
This leads to the two mean field equations 



2 " 

k 

2C5'(5'+l) 



" C"l/3k| , C+|/3k| 
Ak,i 



Ak,2 

1 1 

Ak.i Ak,2 



, (15 



(19) 



Using the values of s and fi thus obtained, we calculate 
the gap to triplet excitations. The dimer-Neel transition 
occurs when the triplon gap vanishes at J± = J±c- Using 
the above equations, we arrive at the following two results 
at the critical point, (i) The value s at the dimer-Neel 
critical point is independent of spin and given by 



2/V, 



■E 



I7k| 



6 



I7k| 



6 



V9-3|7k| V9 + 3|7k| 



(20) 



A numerical evaluation shows Sc ~ 0.872. (ii) We find 
the location of the dimer-Neel critical point 



S{S 



1) N^^ 



1 



1 



v/9-3|7k| v/9 + 3|7k| 



(21) 



A numerical evaluation yields J±c ~ 1.748S'(S' + 1). For 
5 = 1/2, the mean field resuh, J^c[S = 1/2] w 1.311, is 
somewhat smaller than the QMC value. For higher spin, 
the mean field critical points, J±c[S — 1] ^ 3.496 and 
J±c[S = 3/2] w 6.555, are significantly smaller than the 
corresponding QMC results. This is summarized in Ta- 
ble ini Remarkably, as shown in Fig El the scaling result 
J±c ^ 5(5' + 1) from mean field theory appears to be 
valid even for the exact QMC results on the honeycomb 
lattice. We have also explicitly checked that triplon con- 
densation of the mode with energy Ak,i at momentum 
k = (0,0) yields Neel order on the honeycomb bilayer. 
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V. BEYOND MEAN FIELD THEORY: 
VARIATIONAL ANALYSIS 

Corrections to the mean field Hamiltonian arise from 
triplet-triplet interactions, and coupling to higher spin 
objects such as quintets and heptets. As a function of S, 
we find two regimes where two different correction terms 
dominate. For S — 1/2, the only correction stems from 
triplet-triplet interactions since higher spin states are ab- 
sent. For S > 1/2, the dominant correction arises from 
coupling to higher spin (quintet) states. Ordinarily, such 
quintet terms can be ignored as the energy cost of ex- 
citing quintets is large; however, these terms scale as S'^ 
as opposed to the S'^ scaling of the triplet-triplet inter- 
actions and they play an increasingly important role for 
larger S. These two correction terms are separately dis- 
cussed in the following two sections. 

Specifically, for the two regimes S" = 1/2 and S > 1/2, 
we identify the leading correction term and take it into 
account using a variational approach. With the lead- 
ing correction, the Hamiltonian takes the form Hj^^^^ — 

if(o) + A7?^ and ^ H^^^+AHq- We treat AH as a 
perturbation acting upon the states of H'^'^\ As a varia- 
tional ansatz, we assume that the effect of the correction 
terms is entirely accounted for by a renormalization of the 
parameters s and fj, which enter the mean field Hamilto- 
nian, i/^"^ or Hq\ We choose /x to enforce single boson 
occupancy per site on average. The perturbations Ai7, 
for both regimes, preserve total boson number. Thus, it 
suffices to evaluate total boson number using This 
gives us the constraint 



(22) 



where the expectation value is evaluated with respect to 
(For the honeycomb lattice case, there is an addi- 
tional sum over the sublattice degree of freedom in the 
above equation). This leads precisely to the mean field 
number constraint in Eq. [TT] or Eq. I18[ which can now 
be used to determine /Lt. The parameter s is chosen to 
minimize the ground state energy, evaluated to leading 
order in perturbation theory. For S = 1/2, we find that 
the leading correction is obtained within first order per- 
turbation theory in AH. For S > 1/2, the dominant 
perturbing terms require us to go to second order in per- 
turbation theory. In the next two sections, we discuss 
these correction terms in detail. 



only correction. For S > 1/2, this constitutes one term 
in a slew of correction terms. For any spin S, the triplet 
interaction terms are given by 

e {x,y,z} 

We note that there are no cubic terms in triplet opera- 
tors. As described in Ref.'ssi, this makes our bilayer prob- 
lem qualitatively different from other dimerized states 
such as the spin- 1/2 staggered dimer on the square lat- 
tice. Typically, triplet-triplet interactions such as those 
of Eq. [221 are taken into account within a self-consistent 
Hartree-Fock approximationi^i^i Here, we take the in- 
teractions in Eq. [231 to be a perturbation acting on iJ^"'' 
and evaluate the first order correction to ground state 
energy. To this end, we decouple AiJ^ -* using bilinears 
that possess finite expectation values at the level of mean 
field theory: 



(24) 



Here, i and i + 5 are nearest neighbours on the square 
lattice. Explicit expressions for p and A are given in 
Appendix [X] We note that p and A are functions of the 
variational parameters s and The first order energy 
correction due to triplet interactions is given by 



r(t) 



(25) 



Thus, the variational energy of the ground state upon 
including the triplet interaction term is given by 



(26) 



where E^"^ is as defined in Eq. 1101 The parameter s is 
chosen to minimize this energy. We find that the triplet 
interactions reduce the stability of the dimer phase and 
shift J±c to larger values. For S = 1/2, this leads to 
a renormalized transition point J±c ~ 2.58, very close 
to the QMC result. For S > 1/2, the renormalization 
is too weak to account for the discrepancy between the 
earlier mean field result and the QMC data. These triplet 
corrected results for the square lattice are summarized in 
Table n 



B. Honeycomb 



VI. TRIPLET INTERACTION CORRECTIONS 

A. Triplet Interactions on square lattice 

Staying within the singlet-triplet sector, the term we 
have ignored in the mean field treatment is the triplet- 
triplet interaction term. For 5* = 1/2, there are no higher 
spin bosons beyond the singlet-triplet sector, so this is the 



The interaction between triplets on the honeycomb lat- 
tice is given by 



AH 



(t) 
o 



A,w X 



M, U, t/J, u', w' 



e {x,y,z} 
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The operators 5 are such that the sites {i, A) and {i+S, B) 
are nearest neighbours. This interaction term contributes 
to the ground state energy at first order in perturbation 
theory. To evaluate this correction, we quadratically de- 
compose the interaction using the fohowing two bihnears: 

itlAjl+s,BJ ^ '5.,.A, (28) 

with the expectation vahies to be evaluated using the un- 
perturbed Hamiltonian Hq\ Explicit expressions for p 
and A are given in Appendix|B] The first order correction 
to ground state energy is given by 

Ai?« = ^iV^[p2„A2]. (29) 

The parameter s is chosen to minimize the energy 

<L(^-A^) = <+Ai?g^ (30) 

As on the square lattice, we find that the triplet interac- 
tions reduce the stability of the dimer phase and shift J±c 
to larger values. For S ~ 1/2, this leads to a renormal- 
ized transition point J±c ~ 1.59, which is in reasonable 
agreement with the QMC result J±c = 1.645(1). For 
S > 1/2, however, the renormalization is again too weak 
to account for the QMC data. These triplet corrected re- 
sults for the critical point on the honeycomb lattice are 
summarized in Table HIl 

VII. QUINTET CORRECTIONS 

In the previous section, we have seen that triplet cor- 
rection terms lead to a reasonably good agreement with 
QMC results for the dimer-Neel quantum critical point 
for 5 = 1/2. However, they fail to account for the signifi- 
cant discrepancy between QMC and bond operator mean 
field theory for S > 1/2. This leads to us to suspect that 
higher order spin excitations on the dimer bonds must 
be responsible for this difference. Upon including quin- 
tet terms, the spin operators at a site i contain a large 
number of terms as given in Eq. (20) and Eq. (21) of 
Ref. [23 and reproduced in Appendix [C] for convenience. 

Using these spin expressions to rewrite the Hamilto- 
nian in Eq. [Tl we find that correction terms beyond mean 
field theory, including those involving quintet states, may 
be grouped as 

AH ^Dtttt + sRttgiS^) + FttyyiS^) + G,,,,{S"). (31) 

The subscripts indicate the composition of the terms in 
terms of bond operators. The scaling of each term with 
S is indicated in parentheses. For example, Rttq{S^) is 
composed of terms which involve two triplet operators 
and one quintet operator, and the coefficients of these 
terms scale as S^. The term which we have accounted 
for in the previous section is Dtttt , which scales as S*^ and 



contains four triplet operators. Naively, terms involving 
quintets should be less important due to the energy cost 
of exciting quintets. However, we see from the above 
classification of terms that the coefficients of Rttq{S^) 
and Fttqq{S'^) increase rapidly with increasing spin. We 
find that RuqiS"^) is, in fact, the dominant contribution 
for all 5* > 1/2. (For the case of 5* = 1, we have explic- 
itly checked that this term dominates over triplet-triplet 
interactions encoded by Dtttt - see Table |I]) . The term 
FttqqiS"^) is suppressed because it involves two quintet 
operators which act on different sites. In our variational 
scheme, this term will contribute to the ground state en- 
ergy at second order in perturbation theory. However, 
as the quintets are taken to be localized excitations, this 
term will involve intermediate states with two quintet 
excitations. Therefore, it will contribute much less than 
Rttq{S )■ 

In the vicinity of the dimer-Neel transition, we assert 
that RttqiS"^) win remain the dominant correction term 
for any 5' > 1/2 even if higher spin states such as heptets, 
nonets, etc., are included. As the dimer-Neel transition 
occurs via condensation of triplet excitations, it is rea- 
sonable that the dominant corrections come from quin- 
tets which are immediately higher in energy than triplets. 
Heptets, nonets, etc. occur at much higher energies and 
are unlikely to affect the triplet condensation point. To 
argue that this is indeed the case, we first note that the 
Hamiltonian of EqUcan change the spin of a bond by ±1 
at most (this can be seen from the rotation properties of 
a single spin operator acting on a bond eigenstate). For 
example, if we restrict our attention to one particular 
bond, the Hamiltonian connects a triplet state to singlet, 
triplet and quintet states. The matrix element connect- 
ing the triplet to a nonet state (or a state of even higher 
spin) is zero. Similarly, on a given bond, the heptet state 
has non-zero matrix elements only with quintet, heptet 
and nonet states. The resulting terms in the Hamilto- 
nian involving heptets, nonets, etc. will not contribute 
at second order in perturbation theory, but will only ap- 
pear at higher order. As an illustration, upon including 
heptets, the Hamiltonian can have a term of the form 

m1i,nt] m'^j,^' ■ Clearly, this term does not contribute 
to ground state energy at first or second order. In ad- 
dition, at whichever order it contributes, the energy de- 
nominators will involve large heptet excitation energies 
which will further suppress the energy contribution. 

In summary, in the vicinity of the dimer-Neel transition 
for any value of S > 1/2, the leading correction to bond 
operator mean field theory comes from sRttq{S'^). We 
write 

^jj(s>i/2) _ ^jjid) ^ si?tt,(52). (32) 

Note that if we were to use a path integral approach to in- 
tegrate out the quintet excitations at this stage, we would 
be led to an effective triplet-triplet interaction which is 
enhanced by a factor of S'^ compared to the bare triplet- 
triplet term discussed in the previous section (although 
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it would be divided by the quintet gap which scales as 
near the Neel to dimer transition). Here we follow a 
different route, similar in spirit, and treat this term per- 
turbatively assuming the quintet states to be local exci- 
tations. The energy cost of creating a quintet is given by 
Eq. [3l We measure this energy cost from the Lagrange 
multiplier fi, to get 

£g-At = ./±{3-S'(S' + l)}-Ai (33) 
as the energy cost of a quintet excitation. 

A. Quintet corrections on the square lattice 

The terms in sRttq{S^) may be organized as 



i n— — 2,--- 



S 



(34) 



The operator Tf^^j^^ is composed of triplet bilinears. The 
index 5 sums over the four nearest neighbour vectors on 
the square lattice. The explicit form of these operators 
is given in Appendix iDl The operator Rttq{S^) does not 
contribute to ground state energy at first order, as it is 
linear in quintet operators. The energy correction, at 
second order, is given by 



AE, 



(q)_ 



(ok.^„^(T^"l+,,)V)(^l9L^^:+.|o) 



Eo 



(35) 



The index a sums over all excited states of the 
variational Hamiltonian. The only intermediate states 
that contribute are those with a single quintet. In our 
variational formalism, we take the quintets to be local 
excitations. This constrains us to {i = i'), (n = n'). This 
leaves us with 



Eq — Ey 



(36) 



The intermediate states \v) which contribute involve a 
single quintet excitation. Within the triplet sector, at 
zero temperature, the intermediate states can have either 
(a) no triplon quasiparticles, or (b) two triplon quasi- 
particles. The contribution from states with no triplon 
quasiparticles vanishes due to global spin-rotational sym- 
metry of the Hamiltonian. The energy correction from 
two triplon intermediate states is evaluated to obtain the 
energy correction, \e\^\ The complete expression is 
given in Appendix |DJ 

Being second order in Rttq{S'^), the energy correction 
from quintet coupling naively scales as for large S. 
However, the energy denominator involves the energy of 
quintet states which is proportional to J±. Close to the 



s 


QMC 


MFT 


MFT + 


MFT -f 








triplet interactions 


quintet coupling 


0.5 


2.5220(1)" 


2.287 


2.568 




1 


7.150(2) 


6.098 


6.387 


7.32(14) 


1.5 


13.634(3) 


11.434 


11.714 


13.32(1) 



TABLE I: Value of Jxc on the square lattice from different 
methods for different values of S. MFT stands for mean field 
Theory. The QMC data for S=l/2 is from Ref. [TJ. The 
column 'MFT-I- Triplet interactions' gives variational results 
appropriate for S=l/2. The column 'MFT-|-quintet coupling' 
gives variational results appropriate for S > 1/2. 



dimer-Neel transition, at mean field level, J± approxi- 
mately scales as for large S (see Eq. [T3]) . We expect 
perturbative corrections to preserve this scaling of J±c 
with S^. Thus, near the dimer-Neel transition, AeIj*-^ 
scales as S'^/S'^ ^ S'^. The ground state energy to lead- 
ing order in perturbation theory is thus given by 



p(S>l/2)/- _ p(0) 



(37) 



This energy is a function of s and fi. As discussed earlier, 
fi is tuned to enforce single boson occupancy per site. 



while s is chosen to minimize E, 



(S>l/2) 



{s,^J.). 



Having determined s and /i variationally, we can find 
the gap to triplon excitations as a function of Jj_. The 
Dimer-Neel transition is indicated by the vanishing of 
the triplon gap in the variationally obtained state. As 
summarized in Table HJ the renormalized critical points 
obtained in this manner are within 2.5% of the QMC 
results. While the precise quantitative agreement is per- 
haps fortuitous, and will certainly change depending on 
the nature of the approximations made, the important 
problem we have resolved is to show that the large dis- 
crepancy between QMC and simple bond operator mean 
field theory for 5 > 1/2 can be accounted for by virtual 
quintet excitations. 



B. Quintet corrections on the honeycomb lattice 

On the honeycomb lattice, the terms in sRttq{S'^) may 
be written as 



sRttq = sY^ E 



s 



(38) 



The operators Ai^i+s and lii^i-s are triplet bilinears cen- 
tred on nearest neighbour bonds. We give their explicit 
forms in momentum space in Appendix [E] The terms in 
RttqiS"^) contribute to ground state energy only at sec- 
ond order in perturbation theory. The energy correction 
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s 


QMC 


MFT 


MFT + 
triplet interactions 


MFT + 
quintet coupling 


0.5 


1.645(1) 


1.312 


1.588 




1 


4.785(1) 


3.498 


3.774 


4.80(9) 


1.5 


9.194(3) 


6.559 


6.837 


9.58(18) 



TABLE II: Value of J±c on the honeycomb lattice from dif- 
ferent methods for different values of S. MFT stands for mean 
field Theory. 



may be written as 



a^O i.7i 



5' 



a X 



i,i-\-S 



\0)/{Ea-E^} + {A^ B), (39) 



where the index a sums over all excited states of Hq' . 

As the terms in Rttq{S^) involve one quintet operator, 
only intermediate states with a single occupied quintet 
state will contribute. 



(q) 
O 



{Eq - E,} 



{Eo - E,} 



exchange Jj_c scales as S'(S'-|- 1) within, both, bond oper- 
ator mean field theory and QMC simulations. However, 
there is a systematic deviation between bond operator 
mean field theory and QMC, with the deviation itself 
scaling as ^ S"^. Our variational extension of bond oper- 
ator theory to include the dominant triplet and quintet 
excitations successfully captures this systematic devia- 
tion and gives a more precise estimate of Jj_c. 

Bi3Mn40i2(N03) provides an example of a bilayer 
honeycomb antiferromagnet^^ with S = 3/2, where 
strong interlayer exchange couplings '-^ 2 Ji have been in- 
ferred from electronic structure calculations^^. Despite 
this strong bilayer coupling, our study indicates that 
this material would be deep in the Neel ordered phase 
if there are no other frustrating interactions. We are 
thus forced to attribute the observed lack of magnetic 
order in Bi3Mn40i2(N03) to frustration effects arising 
from further neighbor couplings; such further neighbor 
interactions have been shown to drive a variety of new 
phases on the honeycomb lattice.—"— One recent exam- 
ple of a dimer system with S = 1 is the triangular dimer 
materiaii2Ji Ba3Mn2 08. Our approach could be applied 
to understand the triplon spectrum and the effect of quin- 
tet corrections in this material. In particular, our work 
shows that extracting exchange couplings from fitting ex- 
perimental data to bond operator mean field theory may 
not yield precise estimates. In summary, our work pro- 
vides a starting point to think about the physics of high 
spin Heisenberg antiferromagnets in a variety of model 
systems and materials. 



We evaluate these overlaps in momentum space, as de- 
scribed in Appendix |E1 The intermediate state \v) could 
have either (i) no triplon quasiparticles, or (ii) two triplon 
quasiparticles. However, the contribution from states 
with no triplons vanishes due to global spin rotational 
symmetry. The explicit expression for /S.Eq^ is given in 
Appendix [E] 

Thus, the energy of the ground state to leading order 
in quintet coupling, is given by 
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(S>l/2) ^ (0) 
0,var O 



(q) 
o • 



(41) 



We choose s to minimize this energy. The vanishing of 
the triplet gap in the variationally determined state sig- 
nals the dimer-Neel transition. Our results for Jj_c on 
the honeycomb lattice are shown in Table HT] The renor- 
malized critical points for S ~ 1,3/2 are within 5% of 
the QMC value. 



VIII. DISCUSSION 

In this paper, we have studied the Neel to dimer tran- 
sition in Heisenberg antiferromagnets on bilayer square 
and honeycomb lattices for different spin values using 
QMC and bond operator approaches. The critical bilayer 



Appendix A: Square bilayer: bosonic Bogoliubov 
transformation 



The MFT Hamiltonian of Eq. [7] is diagonalized by a 
pseudounitary matrix. 



jj _ ( cosh 0k sinh^k 
I sinh6'k cosh6'k 



(Al) 



Imposing tanh2(?k = ^ 2ek/(^ + 2ek), we get 
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We have defined new quasiparticle operators given by 
V'k,u = C^k0k,u so that 



tk,u \_ ( cosh6'k smh6'k 
*-k« / \ smh^k cosh6'k 



Tk,u 



(A3) 



The T operators are the triplon quasiparticles. The bi- 
hnears defined in Eq. 1241 may be evaluated using the 
elements of U as follows: 



ik.S 



AN 



k,5 

= 7TrE'(2cos/c, + 2cosfc^)^^4^, (A4) 
4iVj^ , Ak 
k 

^-i]^EK/-k>'''1 

k,<5 

= -^Y. ' {'^ COS + 2 COS ky)^^^^. (A5) 



Appendix B: Honeycomb bilayer: bosonic 
Bogoliubov transformation 



The 1? operators are the triplon quasiparticles. Com- 
pared to the square lattice case, the quasiparticle opera- 
tors have an additional index on account of the sublattice 
degree of freedom. We can express our original triplet 
operators as follows: 

tk,A,u = E (C'kJ^'kJ.u + S'kJl^LkJ.n 
t-k,B,u= (Ck,/^-k./,« + 5kj<jJ.(^4) 

The bilinears defined in Eq. [55] can be evaluated as 
2 



3N± 



C-|/3k| , C+|/?k| 



Ak,i 



A 



k,2 



(B5) 



^ = 5]v[E<*k,A,„^-k,B,j7k 



6Nj_ 



Ak,i Ak,2 



(B6) 



The mean field Hamiltonian of Eq. [15] may be diago- 
nalized by the matrix, 



/ 1 1 0\/Ck,i S-ka \ 



Pk = 



V2 



-6k 6k 
11 



Ck,2 Sk^2 

Sk,i ai 



V -6k bk/\ Sk^2 C7k,2/ 



Here, we have defined 6k = /^k/l/^kl- We take the 
other entries to be hyperbolic functions given by Ck.n — 
coshKk.n and Sk,n = sinhKk.n, with n — 1,2. With this 
definition, this matrix Pk satisfies the pseudo-unitarity 
condition PkcPj^ ~ where a ~ Diag{l, 1, —1, —1}. To 
diagonalize the Hamiltonian matrix Afk, we set 



tanh2Kk4 = /3k/(C-/3k); 
tanh2Kk^2 = -/3k/(C4-/?k)- 



(Bl) 



With this choice, the matrix Pk diagonalizes the Hamil- 
tonian, 

PjlMkPk = Diag{Ak,i, Ak,2, Ak,i, Ak,2}. (B2) 

where Ak,i/2 are as defined in the main body. We trans- 
form the triplet operators defined in Eq. [16] into new 
quasiparticle operators using 



tk,A,u \ 
tk.B.u 
''-k,A,u 



Pk 



/ l^k,l,« \ 
1?k.2.« 
^^-k.l.. 

v^ik,2,«y 



(B3) 



Appendix C: Spin operator expressions including 
quintet terms 

Including triplet and quintet operators, the complete 
expression for the spin operators on the two layers of the 
bilayer are^i 



5i=i,2 = (-1) 



, , /(25-l)(25 + 3) r,. t , ^ 

+ (-l)'y ^ 1 ^[(iI,_ift,-2 -<2ki) 



+ 9l,2*a +'?1-i9j-2, (CI) 

with S~^^^ 2 being its Hermitian conjugate, while 



+ (-1) 



(25-l)(25 + 3). /I t , ,J , ) 



5 3 



(C2) 
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Appendix D: Square Bilayer: inclusion of quintets 

The spin operators with the inclusion of quintets are 
given in Eq. 21 of Ref. HJ. Using this reference, we now 
give exphcit expressions for Rttq{S^). In the main text, 

we defined Ruq(S'^) in terms of triplet bilinears Tffj^g- 

Here, we give expressions for in momentum space. 

We use the Fourier transform convention 



t 



i,uG{x,y,z} 



1 v-^ 
k 



(Dl) 



The operator T^i^g is composed of bilinears of 

the form ti^u{U+s,v ± tl+5.v)- Using the Fourier 
transform, this generic bilinear may be written as 

Ek,p i-k+p.n(ik,« ± ilk, Je^'^-'e^P--^^ 
Thus, we may write 



^ k.p 



"?7k, 



(D2) 



where ryk = '^^e 



jk-5 



2(cosfc2: + COS ky) and the coef- 



ficient M = y^ g(s+i)(25^i)(2g+3) ^ rpj^g explicit forms of 
^-k+P,k are: 

-^-k+p,k ~ ^-k+p,a:(^k,a; + ^-k.^) ^ ^-k+p,i/ (^k.j/ +^^k,i/) 
+ *?-k+p,a;(^k,i/ + ^Lk_j,) + i^-k+p/y (ik.x +^^k,2;) 



-k+p,k — *-k+p,2(ik,a:+i_k,2:) + ^-k+p.a;(tk,2; +i ^k^^) 



*^-k+p,2(^k,y + ^'_k^,.) + *i-k+p,y(ik,z+i_k ; 



T 



[0] 
k+p,k 



2 r 



^— k+p,x (^k,a; — k.a:) 



k+p,k 



T 



f-[2] 



k+p,k 



i-k+p,a(ik,y+iLk.y)+2^-k+p,z(^k,2+t'_k,2) 

^^-k+p,z(^k,3;+iLk,j;) ^ i-k+p,a;(ik,z+^Lk.z) 
^i-k+p,z(ik,i,+^Lk,j,) + *^-k+p,a(^k,2+i^k,2) 
^-k+p,a:(^k,x + ^Lk -r) ^ ^-k+p,t^ (^k.j/ +^Lk.i,) 



We have denoted some triplet operators as t and some 
as t. For the purposes of the square lattice, this distinc- 
tion can be ignored. We will use these same expressions 
in the context of the honeycomb lattice also. For the 
honeycomb case, t and t operators will act on different 
sublattices. 

The energy correction due to coupling to quintets is 
given in Eg. 1361 Using the Fourier transformed expression 
in Eq. ID21 we rewrite the energy as 



.S>l/2_M2s2 



(D4) 



ASS 

;v I 

m=-2,--,2 p 

where p is the momentum of the intermediate state. The 
quantity i?p"' is given by 



4"' = E- 



(HEk^l"il+p,k^?k| 

En — E,y 



(D5) 



Here, (— p) is the momentum of the intermediate state 
li/). As described in the Section fVII Al the intermediate 
states that contribute have two triplon quasiparticle 
excitations and one quintet excitation. Within the triplet 
sector, an intermediate state with momentum (— p) may 
be represented as 



W2-triplo7l) — ■^q-p,u''''lq,«'|0)- 



(D6) 



With this parametrization, the sum over intermediate 
states {v) may be written as 



E--E E 



(D7) 



Evaluating the matrix elements using this parametriza- 
tion of the intermediate state, we find that the energy 
contribution is the same from every m-sector, i.e.. 



^ ^— k+p.a; (^k,y 



*i-k+p.j/(^k.x+^lk ~ ™- quantity Ep is given by 



Er, 



^E 



[sinh2(6»q)772_q{cosh(26lp_q) + sinh(26lp_q} + sinh2(6lp_q)?72{cosh(26'q) + sinh(26lq)} 



(D8) 



p+q 



Appendix E: Honeycomb bilayer: inclusion of 
quintets 

In the main text, we defined RttqiS"^) in terms of triplet 
bilinears and 3^}^^- Here, we give expressions 



for v4-"-'_^^ and -Bj"'^, in momentum space. We use the 
Fourier transform convention 

ii,ae{A,B},ue{x,y,z} — , ^g.k.ite''^''"' ■ (El) 



13 



(i) The terms in ^["+5 are of the form ti^AM{ii+s,B.v + where 



\+SB v)- Using our Fourier transform convention, we 
te 

M 



may write 



^-k+P,ke 7k, 



k,p 



where 7k = 
coefficient M - 



ik-<5 



= 1 + e 



+ e 



— ika—ikh 



(E2) 



and the 



S(S+l)(2S-l)(2S+3) 



30 



is the same as that 



defined for the square lattice case. The exphcit forms of 
^[^j^^p ^ are the same of those of T^l._^_p ^ given in Eq 
with the foUowing redefinition: 



^k,ii — iB,ic,u 



(E3) 



(ii) The terms in B^"^_g are of the form ti^B.u{ti-s.A,v^ 

^1-5 All)- Using our Fourier transform convention, we 
write 



^1,1-6 - 



M 



' k,p 



'7-k 



(E4) 



Exphcit expressions for B_l_^p ^ are the same as those 
given in Eq. ID3I but with the foUowing redef- 



of ' 
. -k+p,k 

inition: 



^k,ii — ^S.k.u 
tk,u — ^j4,k,u 



(E5) 



The quintet energy correction on the honeycomb lattice 
may be rewritten as 



N±/2 



(E6) 



IHEki':tp,k7k|0)| 



Eq — E^ 



The only intermediate states {v) that contribute to the 
energy are states with two triplon quasiparticle excita- 
tions. A generic intermediate state with momentum (— p) 
may be characterized as 



W2-triplon) — ^lqj,«^q-p,g.ii|0) 



(E8) 



Using this parametrization of a generic state, the sum 
over intermediate states in Eq. IE7I becomes 



(E9) 



u^O q f,g£{l,2}u,v£{x,y,z} 

There is a factor of 1/2 to account for double counting as 
(q' = p — q, /' = g,g' ~ /) corresponds to the same state 
as (q, /, 5). Evaluating the necessary overlaps, we find 
that the contribution from each m is the same (^[T*') = 
(Sjr') = -Bp for m = -2,--- ,2. The quantity Ep is 
given by 



q,/,9 



[5q,/(-l)f (5-p+q.3 + Cp-q.3)|7p-q| + ^p-q,g ( - 1 ) ^ (^-q,/ + Cq,/)|7q| 



£5 - /i + A_qJ' + Aq_p.g 
I 



(ElO) 



By plugging these expressions into Eq lE6[ the correction to ground state energy may be computed. 
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